function U = poisson_solver_2d(mesh, pde)

A1 = assemble_matrix_A_2d(mesh.node, mesh.elem, pde.cor, mesh.basis_type, "dx", "dx");

A2 = assemble_matrix_A_2d(mesh.node, mesh.elem, pde.cor, mesh.basis_type, "dy", "dy");

A = A1 + A2;

b = assemble_vector_b_2d(mesh.node, mesh.elem, pde.f, mesh.basis_type, "x");

[A, b] = process_boundary_2d(A, b, mesh, pde);

U = A\b;

end